# Merge building footprints to parcel boundaries
pacman::p_load(tidyverse, sf, cowplot, data.table, qs)

## Preliminaries
rm(list = ls())

source("code/globals.R")

fipslist <- qread(file.path(WORKING, "fips_to_merge_to_fprints.qs"))

fipsconcordance <- data.frame(
  statename = c("Arizona", "California", "Colorado", "Oregon", "Washington"),
  statefips = c("04", "06", "08", "41", "53")
)

# For a given state and footprint dataset, load all footprints and subset to the desired counties in that state
run_bldg_subsets <- function(statefipscode, msftrelease) {
  # statefipscode <- "53"
  # msftrelease = "2021release"
  
  msftfilename <- fipsconcordance$statename[fipsconcordance$statefips == statefipscode]
  
  footprints <- st_read(file.path(INPUT_PUBLIC, "Microsoft_Buildings", msftrelease, paste0(msftfilename, ".geojson"))) %>% # slow (15 mins for CA)
    st_transform(crs = 6339)

  fipsinthisstate <- fipslist[which(substr(fipslist, 1, 2) == statefipscode)]
  
  lapply(fipsinthisstate,
    subset_footprints_counties,
    footprintdata = footprints
  )
}

# Load in county shapefile
counties_shp <- st_read(file.path(INPUT_PUBLIC, "CensusTiger", "tl_2019_us_county", "tl_2019_us_county.shp"))

# Save all footprints for a given county from the passed footprint data
subset_footprints_counties <- function(targetfips, footprintdata) {
  # targetfips <- "53047"
  # footprintdata <- footprints
  # END DEBUG
  
  print(targetfips)
  # ptm <- proc.time()
  
  # Lound boundary for this county
  bufferedcountybound <- counties_shp %>%
    filter(GEOID == targetfips) %>%
    st_transform(crs = 6339) %>%
    st_buffer(1000)
  
  # Trim footprints to within county
  buildings_county <- st_join(footprintdata,
                              bufferedcountybound,
                              join = st_intersects,
                              left = FALSE)

  qsave(buildings_county, file.path(WORKING, "buildings", "countysubsets", paste0("b", targetfips, ".qs")))
  
  return(NULL)
}

run_bldg_subsets("41", "2021release")

run_bldg_subsets("04", "2021release")

run_bldg_subsets("08", "2021release")

run_bldg_subsets("53", "2021release")

run_bldg_subsets("06", "2018release")

# MERGE BUILDING FOOTPRINTS TO PARCEL BOUNDARIES ------

### Function for 1 county
spatial_merge_footprints_parcels <- function(targetfips) {
  # targetfips <- "53047" # START DEBUG
  print(targetfips)
  
  # Load the merged sf object containing homes & parcel boundaries inside fire perimeters
  homes_lot_lines <- qread(file.path(WORKING, "lotlines", "merged", paste0("merged", targetfips, ".qs"))) %>%
    st_transform(crs = 6339)

  ## Note there are multiple parcels for each UAPN-incidentid pair, because of mutliple polygons in parcel data in some cases
  ## I am ultimately going to choose the largest 1 structure on any polygon per PARNO, so this is ok.

  homes_lot_lines <- homes_lot_lines %>%
    dplyr::select(FIPS, PARNO, incidentid, ImportParcelID) %>%
    dplyr::mutate(FIPS = sprintf("%05d", FIPS)) # convert FIPS to correct format
  
  # Account for Unknown WKB type 12
  ckClass <- function(x) "MULTISURFACE" %in% class(x)
  mts <- unlist(lapply(homes_lot_lines$geometry, ckClass))
  if(TRUE %in% mts){
    homes_lot_lines$geometry <- st_cast(homes_lot_lines$geometry, "MULTIPOLYGON")    
  }

  # Subset the buildings layer using these home locations
  areaofinterest <- homes_lot_lines %>%
    st_bbox() %>%
    st_as_sfc() %>% 
    st_sf() %>% 
    st_buffer(1000)
  
  buildings_county <- qread(file.path(WORKING, "buildings", "countysubsets", paste0("b", targetfips, ".qs"))) %>% # load the building file clipped to this county's buffer
    st_join(areaofinterest, join = st_intersects, left = FALSE) 

  ## Calculate building area
  buildings_county$bldgarea <- st_area(buildings_county)

  ## Calculate roof centroids; save centroids and bldg polygons in separate vars
  roofinfo <- buildings_county %>%
    st_transform(crs = 4326) %>% # temporarily so can save bldg polygons as WGS84
    mutate(roofpolygon = geometry) %>% # save the building polygon too - used in wall-to-wall neighbor distances.
    st_transform(crs = 6339) %>% # back to 6339 projected system
    st_centroid() %>% # get roof centroids in a projected system (so they'll be correct; st_centroid doesn't work on lat long)
    mutate(
      lon_roof = st_coordinates(st_transform(., crs = 4326))[, 1],
      lat_roof = st_coordinates(st_transform(., crs = 4326))[, 2]
    ) # Save roof center coords in WGS84

  ## - Join parcels and building roof centroid locations
  parc_roofs <- st_join(homes_lot_lines,
    roofinfo,
    join = st_intersects,
    left = TRUE
  )

  ## Identify biggest building on each parcel
  parc_roofs <- parc_roofs %>%
    group_by(PARNO) %>%
    arrange(PARNO, bldgarea, by_group = TRUE) %>%
    mutate(bldgseq = seq_len(n())) %>%
    mutate(biggest = bldgseq == n()) %>%
    mutate(numfootprints = n()) %>%
    dplyr::ungroup()

  ## Create main output: assigned main roof location for each parcel
  parc_locations <- parc_roofs %>%
    dplyr::filter(biggest == TRUE) %>%
    dplyr::select(FIPS, ImportParcelID, lon_roof, lat_roof, numfootprints, roofpolygon)

  # Find parcel centroids and return as alternate measure.
  # Replace NA roof locations with these, as well. These are parcels with no buildings on them (in the MS data).
  parc_centroids <- st_centroid(parc_locations) %>%
    st_transform(crs = 4326) %>%
    st_coordinates() %>% # Gets roof locations as lat/long (WGS84)
    as.data.frame()
  parc_locations$roof_fail <- FALSE
  parc_locations$roof_fail[is.na(parc_locations$lon_roof)] <- TRUE
  parc_locations$lon_roof[is.na(parc_locations$lon_roof)] <- parc_centroids$X[is.na(parc_locations$lon_roof)]
  parc_locations$lat_roof[is.na(parc_locations$lat_roof)] <- parc_centroids$Y[is.na(parc_locations$lat_roof)]
  parc_locations$lon_centroid <- parc_centroids$X
  parc_locations$lat_centroid <- parc_centroids$Y

  # Go back to WGS84 for the parcel boundaries  https://epsg.io/4326
  parc_locations <- st_transform(parc_locations, crs = 4326)

  qsave(parc_locations, file.path(WORKING, paste0("roofpoints", targetfips, ".qs")))
}

lapply(fipslist, spatial_merge_footprints_parcels)

# Assemble and save a single file
roofpoints <- lapply(fipslist, function(x) qread(file.path(WORKING, paste0("roofpoints", x, ".qs")))) %>%
  do.call(rbind, .) %>%
  st_cast()

qsave(roofpoints, file.path(WORKING, "rooftop_locations.qs"))

print("FINISHED FINDING ROOFTOP LOCATIONS")
